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F.E. Harris has been a significant partner in our work on orbital-free density functional 
approximations for use in ab initio molecular dynamics. Here we mention briefly the 
essential progress on single-point functionals since our original paper (2006). Then we 
focus on the advantages and limitations of generalized gradient approximation (GGA) 
non-interacting kinetic-energy functionals. We reconsider the constraints provided by 
near-origin conditions in atomic-like systems and their relationship to regularized versus 
physical external potentials. Then we seek the best empirical GGA for the non-interacting 
KE for a modest-sized set of molecules with a well-defined near-origin behavior of their 
densities. The search is motivated by a desire for insight into GGA limitations and for a 
target for constraint-based development. 

Keywords: orbital-free kinetic energy, cusp condition, Pauli potential, molecular bind¬ 
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1. Introduction 

1.1. Background 

At the Harris Workshop (10-12 Dec. 2014), the second author spoke about recent 
progress by the Univ. Florida orbital-free density functional theory (OFDFT) group of 
which Frank Harris is a member. Substantially all that work is reported in Refs. mm 
and references therein. Earlier work and voluminous references for context are in two 
review articles from our group m- 

Rather than recapitulate the talk and publications, here we provide a particular per¬ 
spective on what has been learnt. The phrase “near-origin” rather than “cusp condition” 
is a clue to the role the external potential plays in enforcing behavior upon generalized 
gradient approximations (GGAs) for the non-interacting kinetic energy. We present some 
new results on near-origin conditions applied to GGAs. These extend work we did with 
Frank Harris [611718] . Then we explore implications of a generic regularization of the usual 
external potential (from a nuclear array) by empirical determination of the most nearly 
optimal GGA for a set of molecular data. That continues the study of binding in sim¬ 
ple molecules by non-self-consistent OFDFT with key ingredients of the methodology 
introduced in our previous publications [6117118] . These ingredients, besides the near-origin 
analysis of the Pauli potential, include (i) the use of a set of nuclear spatial configurations 
for the same molecule; (ii) the use of Gaussian Kohn-Sham molecular densities as input; 
(iii) so-called A E criterion which enforces binding and (iv) the E criterion which enforces 
correct absolute energies. See also the recent work of Borgoo et al. [5] in which the rela¬ 
tionship between binding and the effective homogeneity of approximate non-interacting 
kinetic energy functional is considered. 

Though our research agenda emphasizes functionals for free energy DFT [10J primar¬ 
ily for use in the warm dense matter regime, here we restrict attention to ground state 
OFDFT. There are three reasons. First, ground state OFDFT is a hard challenge (as his¬ 
tory going all the way back to Thomas HU and Fermi [12] demonstrates). That challenge 
is worsened by going to finite-T (one must devise an entropy functional and incorporate 
the intrinsic T-dependence of other functionals). Third, the ground state approximations 
must be reliable and well-founded if there is to be a sensible T=0 K limit for approximate 
free energy functionals. 

1.2. Basics and Notation 

For context and to set notation, the Levy-Lieb version of the Hohenberg-Kohn universal 
functional is 



S[n] =T[n\ + U ee [n ] 


( 1 ) 


with T[n], U ee [ri ], and n( r), the total kinetic energy (KE), total Coulomb energy (Hartree, 
exchange, and correlation), and the electron density at point r respectively (f drn(r) = 
N ei with N e the number of electrons). Assuming that the sum is bounded below, addition 
of an external potential energy E ext [n\ gives the usual DFT variational principle, 


min{£ [n] + E ext [n}} = E 0 [n 0 ] . 

n 


( 2 ) 


Zero subscripts indicate ground state values. 
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The Kohn-Sham [TB] (KS) rearrangement of E invokes an auxiliary non-interacting 
Fermion system with the same density as the physical system. This raises so-called v- 
representability requirements which we assume to be satisfied. The KS system has KE 
and exchange (X) energies T S1 E x , which enable the regrouping of (pD) into 


E[n] 

E c [n\ 

E H [n\ 

E x [n\ 


■\n] + E H [n] + E x [n ] + E c [n] 

’e[n] ~ E H [n} - E x [n} + T[n] - T s [n] 

dr 1 dr 2 n(ri)n(r ? ) 


( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 


E c is the DFT correlation energy. It often is useful to write E xc = E x + E c . V ee is the 
electron-electron Coulomb interaction. The KS Slater determinant $ s [n] is comprised of 
the orbitals from the KS system Euler equations, 


hKs[n\(pi = £Wi , n(r) = ffipfir)] 2 . 

I 


( 7 ) 


Here the fi = 0,1,2 are occupation numbers in the non-spin-polarized case [TTH%] . The 
KS potential is 


vks 

Vh 


Vh 


dr 


Crc T V ex t 
n{ r 2 ) 

27-T 

T “ r 2 


^ err 


5Ex 

5n 


( 8 ) 


For now we leave v ext unspecified. Finally the KS KE is 


N e „ 

T s[Wi}} = J^/i / dr|V0. 
i= 1 7 


(r)| 2 := / <ir(„.j[n(r)] 


( 9 ) 


in Hartree atomic units. This positive-definite integrand form of T s is preferable for 
OFDFT because the integrand of the ordinary Laplacian form of T s can have both signs. 
The two forms differ by a surface integral which is zero for physically significant systems. 


1.3. Essential Challenge of OFDFT 

Posed succinctly, the OFDFT opportunity is that the computational costs of direct 
minimization of Eq. (J3]) scale with system size. In contrast, solution of the KS equations, 
©. has computational cost scaling as IV l or worse. 

The OFDFT challenge can be stated succinctly too. E x [n] is defined in terms of the KS 
orbitals, hence is known exactly only as an implicit functional of n. T s [n] is an implicit 
functional as well. E c [n] is defined in terms of those two. One might try reversion to 
E[n] for construction of approximations but most rigorous knowledge about £ (scaling, 
bounds, limits, etc.) is in terms of the KS rearrangement. Roughly a half century of 
effort has been devoted to Ending good approximations to E x and E xc . Abandoning the 
KS decomposition would discard that resource and, worse, disconnect the result from 
a huge literature of calculations with such functionals. And T s has several rigorously 
demonstrable properties which serve as stringent constraints on approximations [TI5117II5] . 

In short, retention and use of the KS decomposition is practically inescapable. OFDFT 
thus aims at reliable approximations for KS DFT quantities without explicit dependence 
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on the KS orbitals. The allowed variables therefore are n and its spatial derivatives. 
For E xc the consequence is a restriction to the meta-generalized-gradient approximation 
(mGGA) rung of the widely quoted Perdew-Schmidt Jacobs’ ladder of complexity [19]. 
mGGAs depend upon n, | Vrz|, V 2 n, and the KS KE density t orb . Immediately the OFDFT 
challenge is in play: an explicit functional for t orb is required. 

In practice, the highest spatial derivative dependence that so far has been useful for 
t or b is a GGA, to wit 

T GGA [n] = c TF J drn 5/3 (r)ib( S (r)) 

Ctf = lo^ 273 ' ( 10 ) 

F t is called the enhancement factor. For F t = 1, T GGA — Ttf, the Thomas-Fermi func¬ 
tional. The dimensionless reduced density gradient is 

s . = 1 l Vw l - r J Vw l (in 

2(3tT 2 ) 1 /3 U 4 /3 “ n 4/3 ' 1 J 

Remark: The s variable occurs in GGA X functionals also. They have the same form 
as Eq. ([TUP but with n 4 / 3 rather than n 5//3 , F x rather than F t , and a different prefactor, 

Eq. (HDD is a one-point GGA functional. In the OFKE literature there is extensive work 
on two-point functionals, generically f drdr'n a (r)K(r, r')n^(r'). See Section 2.3 of Ref. 
[4] for brief discussion and references. One readily can imagine constructing a GGA for 
two-point functionals but we are unaware of effort along that line. Instead the emphasis 
has been on constructing A'(r, r') via constraints, mostly to match response properties of 
the weakly perturbed homogeneous electron gas. 

Motivated to optimize computational performance, our group has focused on Ending 
and exhausting the limits of single-point GGAs for T s . The remaining discussion assesses 
what we have found, with a focus on the surprising non-universality of approximate OFKE 
functionals, implications for their common use with external potentials of regularized 
Coulomb form, and an empirical attempt to ascertain the limits of GGA performance for 
a particular kind of regularized potential. 


2. Qualitative Distinctions among GGAs for t orb 


Our work began [6] by testing multiple published T GGA functionals. A rough classifi¬ 
cation introduced then was standard GGA and modihed-conjoint GGA (rncGGA). The 
latter term stems from conjoint functionals [2U] . i.e. those for which F t oc F x . Standard 
GGAs include the second-order gradient approximation (SGA) 


rpSGA 

T w [n\ 


Ttf + g T w 


dr 


|Vn(r)p 

n(r) 


( 12 ) 

(13) 


and the von Weizsacker KE T w itself, along with most of the GGAs of the modern era, e.g. 
that by Perdew [21], the PW91 KE functional based on the Perdew-Wang X functional 
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[22]. and those from DePristo and Kress [23], Thakkar [24], and Tran and Wesolowski 
ca, and the APBEK functional based on the PBE X functional [26]. Modified conjoint 
GGAs arise from altering or refining the conjointness conjecture (which is not strictly 
correct; see Ref. 0)- These include Ttf T T\y and our PBE2, KST2 [6]l8] . and VT84F [3] 
functionals. 

The two functional types have qualitatively different performance. Ordinary GGAs 
predict the KE order of magnitude correctly but fail to give binding for simple molecules 
and solids. There are some exceptions for solids in which a pseudo-density is used. mcG- 
GAs do bind simple molecules as well as many solids at least semi-quantitatively, but they 
overestimate the KE strongly. As a consequence the total energy also is strongly overes¬ 
timated (too high). We also have found that these functionals exhibit peculiar sensitivity 
to the type of pseudo-potential used, behavior found by others as well [27] . 

The main difference between ordinary GGAs and the mcGGAs is the enforcement of 
positivity constraints on the mcGGAs. Enforcement is via imposition of requirements 
upon the density in the case that the external potential is Coulombic, 

v ex t(r) = | \ | (14) 

a 1 1 

with Z a the atomic number of the nucleus at site R Q . Such approximate T GGA functionals 
therefore are not guaranteed to be universal, even though T s is. Two questions then arise. 
Are both the overly large KE and sensitivity to pseudo-potentials of mcGGAs connected 
with this non-universality? Is there an example of a T GGA that has both the good KE 
magnitudes of an ordinary GGA and the good binding properties of an mcGGA? We 
address these two issues in the remainder of this paper. 

3. Positivity and Near-origin Conditions 

The Pauli-term decomposition 


T s [n] = T w [n] + T e [n\ 
provides a rigorous bound 
T e [n] > 0 , 


(15) 


(16) 

because 2V is a lower bound to the KS KE [33I34II35I36] . TV also is the exact T s for 
one electron, a fact that will become useful shortly. (It also is exact for a two-electron 
singlet.) The Pauli term potential also is rigorously non-negative: 


, 5 ST e [n] 

W) : = -rTT > 0 
on{r) 


Vr 


(17) 


These are universal properties of T s . For a GGA, the Pauli term separation corresponds 
to a T gga with energy density tg and enhancement factor 


F ,W = F,(s) - . 


(18) 
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Though Tg > 0, it is not necessarily the case that the associated Pauli-term energy 
density tg obeys the same positivity tg > 0 because energy densities are defined only up 
to additive functions which integrate to zero. Refs. [311136M37] chose the canonical form 
for tg (i.e. that which comes from the KS equation), which is positive semi-definite. We 
adopted that argument in Ref. [38]. The consequence, to which we return in Sect. 4, is 

Te(s(r)) > 0 , Vr . (19) 

To have enough additional constraints to determine a useful approximate Fg, we used [3] 
I61I8TI38] requisites of physical many-electron systems, i.e., those with an external potential 
given by Eq. (H. Non-universality enters. 

The nuclear-cusp condition [39] density 

n(r) ~ e~ 2Zr = (1 - 2 Zr) + 0(r 2 ). ( 20 ) 

gives v GGA (r) ~ a/r where a is a constant which depends on the specific enhancement 
factor [ 8 ]. So far as we know, the first mention of this consequence was by Levy and 
Ou-Yang; see the latter part of Section III of Ref. [3Tj 

The one-electron character of the tail region of a many-electron atom [37J forces the 
approximate functional to go over to T\y in that region [TB], For a GGA therefore, 

lim F e (s) = 0 . (21) 

s—>• oo 

Kato cusp behavior Eq. (12U]1 is not exhibited by any density that results from a reg¬ 
ularized potential, e.g., a pseudo-potential. See for example, Eq. ( 6 ) in Ref. [40] and 
associated discussion. Removing that cusp to allow use of compact basis sets (especially 
a plane-wave basis) is the motive for pseudo-potentials. Densities from expansion in a 
finite Gaussian-type basis set, even in all-electron calculations that use Eq. id. also do 
not have Kato cusp behavior. Similarly, the proper tail behavior, also exponential, is not 
found in any finite Gaussian expansion density. Here we focus on the former issue, the 
near-origin behavior of atomic-like systems. 

Consider one-center fV e -electron densities of the flexible form 

n/(r) := H/exp(—Ar 7 ) , l<y <2 

= W7A 3 / 7 

f 4 vrT( 3 / 7 ) ’ 

The norm follows from Ref. [41], with the usual T function. With 7 = 1, A = 2 N e , 
N e = 1, this is the H atom density in the central field approximation. For 7 = 2 it is pure 
Gaussian. For use in what follows, the von Weizsacker potential for densities of this form 

is 

v w = = ^r 7 _ 2 [ 2(7 + 1) - A 7 ?- 7 ] . (24) 

on 8 

With densities of the form (E2D we can explore two simple but illuminating issues. The 
first is to determine the external potential that corresponds to the given density for the 
case N e — 1. Recall the bijectivity of the external potential and the density guaranteed 


( 22 ) 

(23) 







7 


by the first Hohenberg-Kohn theorem [13]. The central-held hydrogenic case is obvious 
but it is instructive to do it in the context of OFKE functionals. The Euler equation is 


5(S + E ext ) 

Sn 


v w + v e + v H + v xc + v ext = ii , 


(25) 


with /i the Lagrangian multiplier for charge normalization. T\y is exact for the one-electron 
case, so vg = 0. Exact exchange cancels the Hartree self-interaction, so vh = —v x , and 
there is no correlation, v c = 0. The von Weizsacker potential (1241) for the hydrogenic 
densities (7 = 1, A = 2N e ) is 


N e N 2 e 

V W —- 7T 

r 2 


For H, fi = — 1, N e = 1, (1251) gives the expected result: 


0 = \-\ + Vext{r)-{-\) 


Vextir) =-• 

r 


Redoing the argument with 7 = 2, N e — 1 gives 


VW = ^(3 - Ar 2 ) 


3A, 


Veit = |A T- +(/!-—), 


(26) 


(27) 


(28) 


the expected quadratic dependence for v ext . 

This elementary exercise illustrates a significant point for approximate functionals. 
Repeat the argument for 7 = 2 but now with the physically important external Coulomb 
potential imposed and with an approximate Tg functional (not necessarily a GGA; for the 
moment the discussion is general). Then the Euler equation becomes 




1 Q \ 

-- - |AV + y + u7 pro » 


(29) 


The only way this can be satisfied is for there to be an incorrect, i.e. non-zero, v ^ pprox f or 
the one-electron case. 

In the case of pseudo-potentials, the argument runs in reverse. Suppose a pseudo¬ 
potential prescription to be used at the so-called one-electron level, i.e., one electron 
outside the core, and suppose it to deliver the form (|28l) . Assume that one can contrive 
a satisfying approximate functional with the property that for N e — 1, the approximate 
functional respects rigorous constraints for the corresponding pseudo-density. Now shift 
to an all-electron pseduo-potential and shrink the core toward the bare Coulomb potential. 
I 11 an arbitrarily small region around the origin, the pseudo-density will remain harmonic 
but the pseudo-potential in almost all space will be essentially Coulombic, leading to the 
kind of mismatch given in Eq. ([29]) . Even at this level (two pseudo-potentials with the 
same regularization procedure but significantly different core radii and populations), there 
is a lack of universality for the approximate Tg. 

For arbitrary 7 dependence, 1 < 7 < 2, the imputed external potential is 


Vext h V\V h T 


—r 7 2 [Ayr 7 


2(7 + 1)], 


(30) 



(with suitably adjusted p of course). Similar mismatch difficulties will occur for all inter¬ 
mediate 7 values, as will the singularities for 7 7 ^ 2 . 

Now consider GGA functionals with arbitrary N e . The GGA Pauli potential is [8j 


7 GA (7 = Con 2 / 3 {| F „( S 2 )-(| S 2 + 2 p 
with higher-order reduced density derivatives 


8Ft , ,(i 4 \ a^F, 

d(s 2 ) +4 V3 7 d(s 2 ) 2 


p := k 


,V 2 n 4 Vn • (VVn) ■ Vn 

n 5/3 ’ Q ■— K n 13/3 

Evaluation with the flexible density (l22j) yields 
s 2 (r) = K 2 A 2 7 2 r 2 ^ 7 _ 1 ) n 7 2 ^ 3 (r) , 


p(r) = K 2 Ayr 7 2 [Ayr 7 — (7 + l)]n^. 2 ^ 3 (r) , 


and 


q{r ) = k 4 A 3 7 3 r 37 4 [A 7 r 7 — (7 — l)]?ij 4 ^ 3 (r) . 

Except for a negative sign, the coefficient of dFg/ds 2 in Eq. (l3Tj) is 

-s 2 + 2p = 2 K 2 A 7 r 7 ~ 2 [-Ayr 7 — (7 + l)]nj 2 ^ 3 (r) . 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


Notice the singularity at the origin for 7 < 2 . Up to a factor of 4, the coefficient of 
d 2 F e /d(s 2 ) 2 in Eq. (1HB is 


-s 4 - q = k 4 A 3 7 3 r 37 4 


yr 7 + (7 - 1 ) 


n 


-4/3 


(37) 


This is non-singular only for 7 > 4/3. 

For small s, one usually enforces gradient expansion behavior on Fg, 

F e = 1 + as 2 (38) 

and only the first derivative term in vg, Eq. (l3TTh is at issue. After a bit of manipulation, 


p(r) = s 2 


s 2 

7 + 1' 

=> -s 2 + 2 p = 2 s 2 

4 y + 1 


Ayr 7 

3 

3 Ayr 7 _ 


The singularity structure in vg then is evident. The general result is 


v e GA [ n f) — con 2 / 3 < - + <m 2 y 2 A 2 r 2(7 1 } n f 2/3 


Lj 1 3 T U ‘ rx / 1L f 

For convenience the two limiting cases are 

,GGAr„_il _ „ „ 2 /3 J 5 2\2 —2/3 


2(7 + 1 ) 

Ayr 7 


(39) 


(40) 


v e [nf, 7 = 1 ] = con / < - + an 2 X 2 n f 


'4 

Ar ~~ 1 


(41) 
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and 


e GA [ n h~t = 2 1 = <W*/ /3 11 +4aK 2 An / 2/ " [3 - Ar 2 ] j . 


(42) 


The takeaway point is that if one sets out to build an approximation constrained to behave 
properly for 7 = 1 (the physical case), the singularity is inevitable and the near-origin 
positivity is determined by the sign of the gradient expansion coefficient a. Our mcGGAs 
are built to have a > 0. However, if the actual density is regularized and has Gaussian 
form near the origin, then if that density is “cuspy” enough, i.e. has large A, even with 
a > 0 and positivity constraints enforced on building the approximation, there still can 
be small-r regions for which v approx < 0 . 


4. Empirical Exploration of the Limits of GGA KE 
4.1. Methodology 

Our approach to the development of GGA OFKE functionals has been to adopt some 
suitable analytical form for the kinetic energy enhancement factor F t with a few pa¬ 
rameters determined from imposing constraints (e.g. correct scaling if applicable, correct 
small-s and large-s behavior) and, if unavoidable, fitting to a small set of training data. 
Interpolation between small- and large -s is defined by the chosen analytical form for F t . 
The analytical forms usually are relatively simple with deliberately limited flexibility to 
avoid introduction of non-physical kinks or other artifacts in that interpolation. In this 
sense, the properties are analogous to those of standard finite basis sets (see for example 
Refs. [2Tl23l2Tl25lfilf7ll3] h 

The unwelcome effects of limited flexibility can be avoided, at least in principle, by use 
of a numerical enhancement factor given on a mesh So = 0, si, ... , s n — s max There is a 
practical barrier however. To determine such a numerical F t requires numerical integration 
in real space of the complicated (n, |Vn|) dependence of the kinetic energy functional 
integrand t or b, Eq. (]9]h evaluated on a numerical s-mesh. Experience demonstrates that 
the result is unphysical, noisy, numerically unstable results. One can see the difficulty 
simply by considering the numerical integration of an s “density of states” on a mesh of 
points sf. 


V( Si ) 



s(r)) « ^2 w jS SiMrj ) , 
3 


(43) 


with Wj the quadrature weights. Numerical experiment shows that a modest change in 
even a very fine r or s mesh (or both) leads to distinctly different results. 

An effective alternative is Pade approximants 02 of high orders such as were used 
recently for analytical representation of common Fermi-Dirac integral combinations [43]. 
They provide the simultaneous flexibility and smoothness required by the numerical in¬ 
tegration in r. Numerical exploration led to the Pade approximant 


Ft(s) 


1 + Ei=l a i st 

1 + EL ’ 


(44) 





10 


of order [9,10] in the variable s (k — 9, l — 10) as a workable compromise between 
flexibility and number of free parameters. 

Only a few parameters in the approximant can be determined from imposition of con¬ 
straints. The remainder must be obtained by fitting. For the present study, the only 
constraint imposed on Eq. (jTTT) is recovery of the correct second-order gradient expansion 
at small s, 



(45) 


This is accomplished by setting a\ — bi, a ,2 — (5/27) + b 2 . To allow maximal freedom 
for the fitted F t we have not imposed the large-s von Weizsacker limit given in Eq. (l2Tj) . 
Studies of X GGA functionals jH] show that the distribution of s, Eq. (J43]), is negligible 
above about s = 3 for most systems of interest, thus suggesting that such large-s behavior 
constraints are not critically important, at least for fitting. Also note well that in what 
follows, we have not imposed any of the positivity constraints, Eqs. (TT6|) . (TTTlh and (lT9]h 
The motivation is to make the empirical fitting as unconstrained as possible. 

The absolute KE versus binding energy dilemma posed at the end of Sect. 2 has im¬ 
plications for the fitting criteria to be used. The usual KE fitting criterion is equivalent 
to the total energy or E criterion [6], namely to minimize the squared energy difference 
between non-self-consistent OFDFT and reference KS energies at system equilibrium ge¬ 
ometries (from standard KS calculations). The procedure is non-self-consistent for the 
OFDFT calculations because KS densities are used as input. The obvious flaw in the E 
criterion (which was investigated in Ref. [8]) is that it forces OFDFT total energies at KS 
equilibrium configurations to be as nearly correct as possible but ignores the shape of the 
KS binding energy curve near equilibrium. 

In Refs. [61IT118] Frank Harris and two of us introduced what was called the A E criterion. 
In it, the objective function to be optimized is formed from energy differences between a 
point away from equilibrium and the equilibrium point as predicted by the reference con¬ 
ventional KS calculation. The objective function has two of those energy differences, one 
from OFDFT, the other from the KS calculations. Obviously the A E criterion enforces 
binding upon the OFDFT approximation while leaving the total energy uncontrolled. The 
result can be an excessively high total energy. 

In the present work we address these two limitations by making a convex sum of average 
versions of the two criteria. The averages are calculated over all atoms, molecules (and 
their geometries) in a training set. To put the two criteria on the same scale, we use the 
mean absolute relative error (MARE) of energy differences rather than average absolute 
energy differences, 



M,i^e 


(46) 


Here, for the nuclear spatial configuration i of molecule iff, A E M}i = E Mti — E M}£ , with 
E M>e the energy associated with the equilibrium nuclear configuration as predicted from 
conventional KS computations, and N the total number of terms in the sum Eq. (l46jh 




11 


Similarly, the mean absolute relative error of the total energy is 



(47) 


The objective function is a convex combination of both 


uj(a) = au E + (1 - a)uAE , 


(48) 


with cl G [0,1]. Minimization of cu(0) is essentially the A E criterion, and conversely 
for w(l). One expects, or at least hopes, that some intermediate a will provide a KE 
functional with both reasonable binding and reasonable absolute energy errors. 

The training set we used includes nine molecules comprised of Erst- and second-row 
atoms and of diverse bonding types along with three closed shell atoms, M = {LiH, CO, N 2 , 
LiF, BF, NaF, SiO, H 4 SiO, H 4 Si 04 , Be, Ne, Ar}. A set of six bond lengths was used for each 
molecule. Molecular geometries were changed by varying the single bond length in the 
diatomics, the central bond length R( Si-O) in H 4 SiO, and by varying A(Si-Oj) in H 4 Si0 4 
deformed in the mode. This set is small by comparison with the training sets used in 
the Minnesota series of XC functionals [45] because our purpose is different. We do not 
seek a broadly useful empirical functional. The issue here is narrower, namely whether 
there exists an OFKE functional which does well both on absolute energies and binding 
even on a small sample of systems. 

One other technical point is that the enhancement factor F t which results from fitting 
is checked for poles on the interval s G [0,1000]. If the denominator of Eq. (1441) has a root 
on that interval, the corresponding set of parameters is rejected. 

All reference KS calculations were done in the local density approximation (LDA) for 
XC (see Refs. [46114714811491150851152153] ) using a triplc-zeta Gaussian-type basis with polar¬ 
ization functions (TZVP) [5411551156] . Orbital-free kinetic energy integrals were calculated 
by numerical quadrature, as in our previous work [6]. Weight functions, w/(r), local¬ 
ized near each center with the properties that w/(r) > 0 and ^ 7 w/(r) = 1 are used to 
represent the multicenter integrals exactly as a sum of atom-centered contributions EH 



(49) 


7=1 


Radial integration of the resulting single-center forms was via a Gauss-Legendre proce¬ 
dure, while integration over the angular variables used high-order quadrature formulae 


[SB]- A dense mesh consisting of 150 radial and 434 angular grid points was used to cal¬ 
culate atom-centered integrals. These computations used routines developed by Salvador 


and Mayer [59] and included in their code fuzzy. 

Before proceeding to results, one should note the implications of these numerical pro¬ 
cedures. The finite Gaussian-type basis inexorably yields Gaussian near-origin behavior 
of the density. Yet the calculations are all-electron in the bare Coulomb external poten¬ 
tial, dm) . This is precisely the inconsistency between external potential and near-origin 
density behavior discussed in Sect. 3. 
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Figure 1. Minimum co(a) and corresponding uae and uje values as functions of a. 


4.2. Results 

There are seventeen independent parameters left in Eq. (144|) after constraining to the 
second-order gradient expansion. Those were optimized to minimize the objective function 
u>(a). Figure |T] shows the u>e and ujae MAREs, Eqs. (1471) and (l46|) respectively, as 
functions of a. The minimum uj(a) value also is shown. It decreases monotonically from 
41% to 0.12%. Up to about a = 0.97, u>e decreases slowly with a few jumps (from about 
5% to 0.8%), while ujae is almost flat (from 41% to 44%). Unsurprisingly, they diverge 
as a —> 1 (ujae = 140%, co = 0.12%), an illustration of the absolute energy vs. binding 
energy dilemma. 

Figure [2] shows the fitted F t and F g for selected a values. Notice the violation of Fg > 0, 
a consequence of the unconstrained fitting. Notice also the structure in F t , especially for 
large a values. 

The right-hand panel of Fig. [2] clearly shows the separation of enhancement factors 
into two groups corresponding to a < 0.99 and a = 1.0. The a = 1.0 curve is almost a 
line (Fg is shown as a function of s 2 ) and is practically indistinguishable from the SGA 
curve (shown for comparison) for s < 1.5. Fg (and F t ) has some structure (oscillations 
around the SGA curve) for s < 1.5 and a < 0.99, though that may change upon changing 
the training set and/or the analytical form for the enhancement factor. Also, during the 
optimization process we noted the existence of many significantly different enhancement 
factors which cannot be discriminated clearly by the objective function uj(a). 

Any choice of 0.1 < a < 0.98 corresponds to the required KE functional, namely one 
which provides semi-quantitative binding, uae ~ 44%, and reasonable absolute energy 
error, ( oje ~ 1%). Table [Tj lists the uje and ojae values for the systems from the training 
set corresponding to the optimized u>(a = 0.95). The highest ujae = 85% corresponds to 
LiH. Values for the optimized uj(a = 1.0) are also shown for comparison. All ujae values 
for a = 1.0 (except the LiH molecule) are between 110% and 200%, signifying no binding. 

The quality of energy curves for the CO and H 4 SiO molecules corresponding to the 
a = 0.95 functional is shown in Fig. [3J The a = 1 curves (functional parameters fitted 
to optimize u(a = 1), i.e. pure E-criterion) are shown for comparison. They have no 
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Tabic 1 

MARE values uje and uj^e (in %) calculated individually for each system, with parameters 
from minimization of u(a = 0.95) and u{a = 1.0). 


System 

a 

= 0.95 

a = 

1.0 

UJE 

wa e 

UlE 

wa e 

LiH 

4.2 

85 

0.2 

85 

CO 

0.2 

28 

0.3 

150 

n 2 

0.2 

35 

0.3 

150 

LiF 

0.1 

80 

0.1 

200 

BF 

0.1 

15 

0.1 

160 

NaF 

0.1 

25 

0.03 

180 

SiO 

0.02 

32 

0.04 

120 

H 4 Si0 4 

0.3 

41 

0.04 

120 

H 4 SiO 

0.02 

50 

0.6 

110 

Be 

4.0 

- 

0.2 

- 

Ne 

0.1 

- 

0.02 

- 

Ar 

0.7 

- 

0.1 

- 

average 

0.83 

44 

0.12 

140 


minima. For the CO molecule, the minimum at a = 0.95 curve is too shallow compared 
to the reference KS result. In contrast, for H 4 SiO the agreement between the a = 0.95 
optimized functional and KS results is excellent. Fig. ® shows a similar comparison for 
two molecules not in the training set. For the simplest, H 2 , both functionals (a = 0.95 
and a = 1.0) provide qualitatively correct binding, but the minimum is too shallow and 
there are large discrepancies at the dissociation limit relative to the KS result. Single- 
bond stretching in H 2 0 is described qualitatively roughly correctly only by the a = 0.95 
functional. As expected, the a = 1.0 functional fails to give binding. 

Table [2] lists the set of coefficients in the kinetic energy enhancement factor Eq. (|44|) 
for the optimized co(a = 0.95). 

5. Summary Discussion 

A GGA OFKE form is the simplest one-point functional which explicitly includes ef¬ 
fects of electron density inhomogeneity. The GGA form, with parameters determined by 
constraints, has been very successful for approximate XC functionals, to the point that 
such functionals dominate in practical calculations. A crucial distinction with respect to 
OFKE is the order of magnitude of the XC energy, about 10% of the total energy. The 
ground state KE, however, has the same order of magnitude as the total energy and the 
KS KE is a large fraction of the total KE. Hence a relative error of even a few percent in 
an approximate OFKE functional will have much bigger impact on calculated properties 
than would the same relative error in an approximate XC functional. 
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S 


Figure 2. Non-interacting kinetic energy (left) and Pauli term (right) enhancement 
factors, F t and F g , as functions of s and s 2 respectively. 



R(C-O) (A) 



Figure 3. Total energy as a function of bond distance for the CO (left) and H 4 Si0 4 
(right) molecules from the training set obtained from a KS LDA calculation and from the 
post-KS orbital-free calculation with approximate GGA functionals. 
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Figure 4. Total energy as a function of bond distance for the H 2 (left) and H 2 0 (right) 
molecules (neither in training set) obtained from a KS LDA calculation and from the 
post-KS orbital-free calculation with approximate GGA functionals. 


This distinction has significant implications for the constraint-based development of 
GGA OFKE functionals. There is a related, but perhaps less-obvious distinction. Because 
exchange in physical systems is purely Coulombic and correlation (as defined in DFT; 
recall Eq. (B|) is mostly Coulombic (the KE contribution is small), it is eminently sensible 
to impose Coulombic constraints on an approximate E xc . The resulting functional should 
be applicable to a broad range of v extj if not truly universal. Experience shows that to be 
the case. Good GGA E xc functionals deliver variations in MARE over classes of properties 
and types of bonding but they are broadly applicable. 

Generating useful constraints on a GGA OFKE functional that preserve universality 
has not been as straightforward so far as for the XC functionals. One way to see the 
underlying difficulty is that T s is a non-interacting system quantity, whereas E xc is, as 
just remarked, predominantly a Coulombic quantity. This means a lack of specificity 
for T s compared to E xc . Surmounting that lack is what we have done by abandoning 
universality and imposing conditions that follow from Coulombic v ext . Precisely because 
they are non-universal conditions, use of pseudo-potentials or Gaussian-type basis sets 
immediately introduces inconsistencies. In Section 3 we have shown how the consequences 
can be delineated clearly with simple one-center densities. 

Section 4 then considered whether any GGA can have acceptable errors in both total 
energies and binding energies. The result is not as encouraging as we would like. For 
a small training set, the best empirical OFKE GGA we have been able to develop so 
far provides a relatively small MARE for the total energy ( oje ~ 0.8%) but only semi- 
quantitative binding (coae ~ 44%). The use of mixed E-AE criteria is essential to get 
both correct total energies and roughly reasonable binding simultaneously. 

The correct second-order gradient expansion was the only constraint imposed on the 
KE enhancement factor. The result is violation of a constraint, Fq > 0, Eq. (fIU]l which 
depends on a particular choice of to- We have not checked whether vq > 0 is violated for 
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Table 2 

Coefficients in Eq. (ITTl) from optimizing co(a = 0.95). 


coefficient 

value 

CL\ — 

bi 

12.100994770272 

Oj2 = 

5/27 + b 2 

10.829496969896 

Q.3 


-27.327919841144 

CI4 


73.841590552393 

a 5 


25.096089580269 

a q 


-45.306369888376 

a 7 


-77.901835391837 

ag 


-20.862438996553 

ag 


67.083330246208 

bi 


12.100994770272 

b 2 


10.644311784711 

bs 


14.896876304511 

b 4 


5.5830951758904 

b§ 


-24.558524755221 

K 


-31.914940553009 

b-[ 


4.3293607988211 

bs 


17.169012815532 

bg 


2.4210601059537 

b\g 


3.2527234245842 


the empirical GGA, but are certain that it will be violated because of the negative slope of 
pempmcai f or g 2 < ^ Q ne might hope that incorporation of more constraints should make 
the functional better or, at least, more nearly universal in the limited sense of improving 
its transferability to different systems and/or conditions. The counterargument is that, 
except for the high-order Pade form itself, the empirical functional was not restricted in 
any other way, a fact which should facilitate optimization (even at the cost of realism or 
transferability). 

A weak caveat is that we have not yet tested the empirical functionals in SCF calcula¬ 
tions. This may be important, because the large -s behavior of the empirical functionals 
for Pade approximants of different orders might be very different, while the Ue and ujae 
errors are very similar. That difference in large-s behavior might be important for SCF 
calculations, but not for post-KS calculations. However, the somewhat disappointing 
performance on post-KS binding energy curves makes this, in our judgment, a somewhat 
problematic conjecture. 

During minimization of u>(a) with respect to the independent parameters of the Pade 
enhancement factor (at fixed a of course), we encountered three main difficulties. First 
is the precision of numerical integration required (very dense radial and angular meshes) 
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to handle roughness in the kinetic energy enhancement factor. Second is that, for any 
given a, the optimization process very frequently sticks in local minima. Optimizations 
therefore were run for multiple values of a. Those results were analyzed to find one or 
a few superior parameter sets (based on values of cje and coae)- Then those parameter 
sets were used as initial ones to start a repeat optimization for all a. Eventually this “by 
hand” procedure yielded optimal sets of parameters for every a. Thirdly, the existence of 
multiple enhancement factors which deliver very similar values of the objective function 
co(a) despite being very different functions of s makes the final functional for each a 
sensitive to the choice of the training set. Increasing the size of the training set might 
help to overcome that difficulty. But we note again that our purpose here is not to attempt 
a general, widely applicable empirical GGA for OFKE but simply to find the best for a 
modest selection of molecules. Even with that narrow goal, the outcome seems to be that 
there are significant limits on what can be expected of a GGA OFKE functional. 
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